function [S ,CPx, CPy, CPz] =  NURBS_surface(CP)
    
    global knot_vector_U
    global knot_vector_V
    global p
    global q
    p=3;
    q=3;
    knot_vector_U=[0 0 0 0.4 0.6  1  1  1];
    knot_vector_V=[0 0 0  0.2  0.5  0.8  1  1  1];
    % CP=[[0,0,0] [0,1,0] [0,2,0] [0,3,0] [0,4,0]   [0,5,0]
    %     [1,0,0] [1,1,1] [1,2,1] [1,3,1] [1,4,0.2] [1,5,0]
    %     [2,0,0] [2,1,1] [2,2,2] [2,3,1] [2,4,0.5] [2,5,0]
    %     [3,0,0] [3,1,1] [3,2,1] [3,3,1] [3,4 0.2] [3,5,0]
    %     [4,0,0] [4,1,0] [4,2,0] [4,3,0] [4,4,0]   [4,5,0]
    %     ];

    %  CPx=[0 0 0 0 0 0
    %       1 1 1 1 1 1   
    %       2 2 2 2 2 2
    %       3 3 3 3 3 3
    %       4 4 4 4 4 4]
    %  CPy=[0 1 2 3 4 5
    %       0 1 2 3 4 5
    %       0 1 2 3 4 5
    %       0 1 2 3 4 5
    %       0 1 2 3 4 5]
    %   CPz=[0 0 0 0 0    0
    %        0 1 1 1 0.2  0
    %        0 1 2 1 0.5  0 
    %        0 1 1 1 0.2  0
    %        0 0 0 0 0    0]
    [m,n]=size(CP);
    for i=1:m
        for j=1:n/3
    %        if j==1
    %         cpx(i,j)=c(i,j)
    %         cpy(i,j)=c(i,j+1)
    %         cpz(i,j)=c(i,j+2)
              CPx(i,j)=CP(i,3*(j-1)+1);
              CPy(i,j)=CP(i,3*(j-1)+2);
              CPz(i,j)=CP(i,3*(j-1)+3);
          % end
        end
    end

     W=[ 1     1     1     1     1     1
         1     1     1     1      1     1
         1     1     2     1     1     1
         1     1     1     1     1     1
         1     1     1     1     1     1];
     [m,n]=size(CPx);
     for ii=1:m
         for jj=1:n
         Pw_x(ii,jj)=CPx(ii,jj)*W(ii,jj);
         Pw_y(ii,jj)=CPy(ii,jj)*W(ii,jj);
         Pw_z(ii,jj)=CPz(ii,jj)*W(ii,jj);
         end
     end

    tag=1;
    for u=0:0.01:1
       for v=0:0.01:1
            tag1=1;
            i=findspan(u,knot_vector_U); j=findspan(v,knot_vector_V);
            for  temp1=i-p+1:i
              N_kpu(tag1)=blending_function_u(temp1,p,u);
              tag1=tag1+1;
            end

            tag2=1;
            for temp2=j-q+1:j
               N_lqv(tag2)=blending_function_v(temp2,q,v); 
                tag2=tag2+1;
            end
            % S(tag,1)=[blending_function(1,3,u),   blending_function(2,3,u),  blending_function(3,3,u)]*[]
            cpx=compute_cp(i,j,p,q,Pw_x);
            cpy=compute_cp(i,j,p,q,Pw_y);
            cpz=compute_cp(i,j,p,q,Pw_z);
            cpw=compute_cp(i,j,p,q,W);
            Sw(tag,1)= N_kpu*cpx*(N_lqv.');
            Sw(tag,2)= N_kpu*cpy*(N_lqv.');
            Sw(tag,3)= N_kpu*cpz*(N_lqv.');
            Sw(tag,4)= N_kpu*cpw*(N_lqv.');
            S(tag,1)=Sw(tag,1)/Sw(tag,4);
            S(tag,2)=Sw(tag,2)/Sw(tag,4);
            S(tag,3)=Sw(tag,3)/Sw(tag,4);
            tag=tag+1;
       end
    end
end
